function vmd_onefluid(fbase, fnum)
%
%   vmd_images(fbase, fnum)
%
%   Reads file with density distribution and file with x-velocity.
%   Creates bunch of images.
%
%   Input:
%        fbase   base of file number (runname)
%        fnum    file number
%
%   Domain information and colormap limits are edited in the file.
%%

%---------------------------------------------------------

% system size

%jx=48; jy=1;  jz=16; jv=256; jdir=2;
jx=48; jy=48; jz=16; jv=256; jdir=8;

% physics paramters

u0=1.24; vPhase=3.53235; uAngle=0;



%-- backround velocities --------------------------------

  theta = 2*pi*(0:jdir-1)/jdir + pi*uAngle;
  ux0 = u0 * cos(theta);
  uy0 = u0 * sin(theta);
  v0  = -vPhase;


  fn = num2str(fnum,'%04d');

  fluid = '.0';

  imgbase = 't100'

%--------------------------------------------------------



  %-- density

  fname = [fbase, '.dfu.', fn, fluid];

  fid = fopen(fname, 'rb');
  f = fread(fid, jv*jx*jy*jz, 'double');
  fclose(fid);

  f = reshape(f, jv,jz,jy,jx);

  dn = sum(f, 1);

  dn = reshape(dn, jz,jy,jx);



  %-- velocities

  fname = [fbase, '.ux.', fn, fluid];
  fid = fopen(fname, 'rb');
  f = fread(fid, jx*jy*jz, 'double');
  fclose(fid);
  ux = reshape(f, jz,jy,jx);
 

  fname = [fbase, '.uy.', fn, fluid];
  fid = fopen(fname, 'rb');
  f = fread(fid, jx*jy*jz, 'double');
  fclose(fid);
  uy = reshape(f, jz,jy,jx);

  fname = [fbase, '.uz.', fn, fluid];
  fid = fopen(fname, 'rb');
  f = fread(fid, jx*jy*jz, 'double');
  fclose(fid);
  uz = reshape(f, jz,jy,jx);



  %-- divergence


  uxx = circshift(ux, [0,0,-1]) - circshift(ux, [0,0,1]);
  uyy = circshift(uy, [0,-1,0]) - circshift(uy, [0,1,0]);

  uxz = circshift(ux, [-1,0,0]) - circshift(ux, [1,0,0]);
  uyz = circshift(uy, [-1,0,0]) - circshift(uy, [1,0,0]);


  uzx = circshift(uz, [0,0,-1]) - circshift(uz, [0,0,1]);
  uzy = circshift(uz, [0,-1,0]) - circshift(uz, [0,1,0]);


  divu = uxx + uyy;


  s = -(uzx.*uxz + uyz.*uzy);


  %-- images ------------

  iz=jz/4+1; iy=jy/2+1;


% colormap bounds

  maxDn    = 0.01;
  maxUx    = max(max(max(abs(ux-ux0(1)))))
  maxUxz   = max(max(max(abs(uxz))))
  maxDivU  = max(max(max(abs(divu))))*0.5
  maxDivUs = max(max(max(abs(s))))


  f = reshape(dn(1,:,:), jy, jx);
  img_slices(f, [imgbase, '.dn.z0'], -maxDn, maxDn)

  f = reshape(dn(iz,:,:), jy, jx);
  img_slices(f, [imgbase, '.dn.z1'], -maxDn, maxDn)

  f = reshape(dn(:,iy,:), jz, jx);
  img_slices(f, [imgbase, '.dn.y0'], -maxDn, maxDn)




  f = reshape(ux(1,:,:)-ux0(1), jy, jx);
  img_slices(f, [imgbase, '.Ux.z0'], -maxUx, maxUx)

  f = reshape(ux(iz,:,:)-ux0(1), jy, jx);
  img_slices(f, [imgbase, '.Ux.z1'], -maxUx, maxUx)

  f = reshape(ux(:,iy,:)-ux0(1), jz, jx);
  img_slices(f, [imgbase, '.Ux.y0'], -maxUx, maxUx)



  f = reshape(uxz(1,:,:), jy, jx);
  img_slices(f, [imgbase, '.Ux_z.z0'], -maxUxz, maxUxz)

  f = reshape(uxz(iz,:,:), jy, jx);
  img_slices(f, [imgbase, '.Ux_z.z1'], -maxUxz, maxUxz)

  f = reshape(uxz(:,iy,:), jz, jx);
  img_slices(f, [imgbase, '.Ux_z.y0'], -maxUxz, maxUxz)


  f = reshape(divu(1,:,:), jy, jx);
  img_slices(f, [imgbase, '.divU.z0'], -maxDivU, maxDivU)

  f = reshape(divu(iz,:,:), jy, jx);
  img_slices(f, [imgbase, '.divU.z1'], -maxDivU, maxDivU)

  f = reshape(divu(:,iy,:), jz, jx);
  img_slices(f,  [imgbase, '.divU.y0'], -maxDivU, maxDivU)


  f = reshape(s(1,:,:), jy, jx);
  img_slices(f, [imgbase, '.divUsource.z0'], -maxDivUs, maxDivUs)

  f = reshape(s(iz,:,:), jy, jx);
  img_slices(f, [imgbase, '.divUsource.z1'], -maxDivUs, maxDivUs)

  f = reshape(s(:,iy,:), jz, jx);
  img_slices(f, [imgbase, '.divUsource.y0'], -maxDivUs, maxDivUs)



return


end

%---------------------------------------------------------
